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The exact ground state of a strongly interacting quantum 
many-body system can be obtained by evolving a trial state 
with finite overlap with the ground state to infinite imagi- 
nary time. In many cases, since the convergence is expo- 
nential, the system converges essentially to the exact ground 
state in a relatively short time. Thus a short time evolved 
wave function can be an excellent approximation to the exact 
ground state. Such a short time evolved wave function can be 
obtained by factorizing, or splitting, the evolution operator 
to high order. However, for the imaginary time Schrodinger 
equation, which contains an irreversible diffusion kernel, all 
coefficients, or time steps, must be positive. (Negative time 
steps would require evolving the diffusion process backward 
in time, which is impossible.) Heretofore, only second order 
factorization schemes can have all positive coefficients, but 
without further iterations, these cannot be used to evolve the 
system long enough to be close to the exact ground state. In 
this work, we use a newly discovered fourth order positive fac- 
torization scheme which requires knowing both the potential 
and its gradient. We show that the resulting fourth order wave 
function alone, without further iterations, gives an excellent 
description of strongly interacting quantum systems such as 
liquid 4 He, comparable to the best variational results in the 
literature. This suggests that such a fourth order wave func- 
tion can be used to study the ground state of diverse quantum 
many-body systems, including Bose-Einstein condensates and 
Fermi systems. 
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We consider a quantum system of N particles with 
mass m described by the Hamiltonian: 
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N 



H = T + V ; T=-\J2^ ; V = 5>(r«) , (1) 

i—l i>j 

where T is the kinetic energy operator, V is a sum of 
pair-wise potentials v(rij) and A = Ti 2 /(2m). 

In imaginary time r = it the many-body time- 
dependent Schrodinger equation can be written as 

-A|*(r)} = ^(r)) , (2) 



with formal solution 
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|d>) = |*(0)> 
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In the coordinate representation ^(R, r) = (R\^(t)) — 
(R\c~ tH \&), where R = {ri . . . rjy} denotes the set of all 



particle coordinates. If the initial wave function |3>) is 
expanded in the set of exact eigenfunctions {<J> n } of the 
Hamiltonian H, then Eq.(3) has the more explicit form, 



*(i?,r) = e 
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Assuming the non-degeneracy of the ground state (E n — 
Eq > ;n / 0), the above wave function becomes pro- 
portional to the exact ground state wave function in the 
limit of infinite imaginary time. A basic strategy is then 
to start with a good trial wave function and evolve it 
in imaginary time long enough to damp out all but the 
exact ground state wave function. 

Since the imaginary time evolution cannot be done ex- 
actly, one usually develops a short-time propagator by de- 
composing e~ rH — e~ T ( T+v * > into exactly solvable parts, 
and further iterate this short time propagator to longer 
time. This is essentially the approach of the Diffusion 
Monte Carlo (DMC) method [1-3]. The need for itera- 
tions introduces the complication of branching, which is 
the hallmark of Diffusion and Green's Function Monte 
Carlo methods [4]. Our idea is to develop a short time 
propagator via higher order decomposition that can be 
applied for a sufficient long time to project out an excel- 
lent approximation to the ground state without iteration. 

First and second order factorization schemes such as 
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are well known, but without iterations, they cannot be 
applied at sufficiently large value of r to get near to the 
ground state. It is also well known that in the context of 
symplectic integrators, the short time evolution operator 
can be factorized to arbitrarily high order in the form 
[5-12] 



-r(T+V) 



n 



e~ a * rT e 



■biTV 



(6) 



with coefficients {a^, b{\ determined by the required order 
of accuracy. However, as first proved by Sheng [13] and 
later by Suzuki [14] (using a more geometric argument), 
beyond second order, any factorization of the form (6) 



1 



must contain some negative coefficients in the set {a,, bi}. 
Goldman and Kaper [15] later proved that any factoriza- 
tion of the form (6) must contain at least one negative 
coefficient for both operators. This means that for de- 
compositions of the form (6), one must evolve the sys- 
tem backward in time for some intermediate time steps. 
This is of little consequence for classical dynamics, or 
real time quantum dynamics, both of which are time- 
reversible. For the imaginary time Schrodinger equation, 
whose kinetic energy operator is the time-irreversible 
diffusion kernel, this is detrimental. This is because 
e -a s rT x e -(r'-r) /(2a,r) j s diffusion Green's func- 
tion. For positve aj, this kernel can be simulated by 
Gaussian random walks. If dj were negative, the kernel 
would be unbound and unnormalizable, with no proba- 
bilistic based (Monte Carlo) simulations possible. This is 
just a mathematical restatement of the physical fact that 
diffusion is an irreversible process. Positive decomposi- 
tion coefficients arc therefore absolutely essential for solv- 
ing any evolution equation having an irreversible compo- 
nent, such as the imaginary time Schrodinger equation. 

Since both classical and quantum dynamics are time- 
reversible, there is a lack of impetus to search for higher 
order factorization schemes with only positive coeffi- 
cients. While higher order factorizations of the form 
(6), with negative coefficients, have been studied exten- 
sively in the literature [10-12], it was only recently that 
Suzuki [16] and Chin [17] have found some fourth order 
(but no higher order) all forward time step decomposition 
schemes. In order to bypass Sheng and Suzuki's proof, 
one must introduce a higher order commutator [V, [T, V] 
in additional to operators T and V used in (6). In this 
work, we use the fourth order factorization scheme [16,17] 
referred to as scheme A: 



G(R,S,t) = (R\e-™\S). 



(10) 



e -r(T+v) = e - 



tT 



-• TV +0(r 5 ) 



(7) 



with V given by 
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V = V+ T -[VAT,V]]=V+ T -2\Y,\V l V\ 2 . (8) 

i=l 



This scheme was also found by Koseloff [18], but his coef- 
ficient for the double commutator term is incorrect by a 
factor of three too large. For a more detailed discussion 
of positive factorization schemes and forward symplector 
integrators, see Ref. [19]. 

To go from state vectors to coordinate wave func- 
tions, we insert complete sets of coordinate states, 1 = 
/ dS \S)(S\ where S = {si . . . s^} and write, for exam- 
ple, the operator equation (3) in the form 

V(R,t) = J dS G(R,S,t) <P(S) , (9) 

where the Green's function G(R, S, r) is given by 



The intermediate coordinates S are also sometime re- 
ferred to as "shadow" positions. Each decomposition 
scheme then corresponds to a specific wave function for 
the ground state. For instance, the first-order scheme 
gives the linear wave function, 

*(R, r) = JdS e~ c < fl - s ) 2 e~ TV ^ , (11) 

where (R — S) 2 = J2iLi( r i — s 2 ; an d where we have 
used the fundamental result that the kinetic evolution 
operator is just the the diffusion Green's function, 

(i?|e^ T |5) oce- c ( fl - s ) 2 ; C = _L . (12) 

4rA 

Similarly, the second order scheme gives the following 
quadratic wave function 

9(R,t) = e-% y W J dS e - c ^ 2 e~i v ^ $(S) . 

(13) 

Finally, the fourth order scheme A produces the following 
quartic many-body wave function, 

*(R,t) = e"^( fl ) J dS'c- c '( R - s ' fe-^(s') x 

J dSc- c '^'- s ^e-^ s ^(S) , (14) 

now with C = 1/(2tA). 

In all these wave functions, there is only a single para- 
mater, the imaginary time r, that we can vary. All else 
are fixed by factorization schemes. If the factorization 
scheme can accurately reproduce the imaginary time evo- 
lution of the wave function, the resulting energy must 
fall monotonically from the initial energy toward the ex- 
act ground state energy with increasing r. To the extent 
that these wave function are not the exact imginary time 
wave function, the energy will eventually rise again. Thus 
for each wave function there is a optimal r where it will 
be "closest" to the exact ground state. 

To test the quality of the above wave functions we use 
them to describe the ground state of a strongly corre- 
lated quantum system of N 4 He atoms interacting via 
a two-body Aziz HFDHE2 potential. [20] At equilib- 
rium, the system is in a liquid state and has a density 
of per 3 = 0.365 (er = 2.556A). The simplest description 
of the ground state is McMillan's Jastrow wave function, 



$(i?) = cxp 
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; u(r) = \(?) . (15, 



with b = 1.2 a. We will us this wave function as our 
initial wave function in all our simulations. 
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For all three wave functions, the expectation value of 
the Hamiltonian can be computed from 



E = 



j dR V(R,t)HV(R,t) 
J dR \$>(R,t)\ 2 ' 



(16) 



The iterated wave functions simply require more integra- 
tion variables. For example, in the case of the linear and 
quadratic wave function, the above can be expressed as 



E 



dR dS L dS R p{R, S L ,S R ) E L (R, S L , S R ) , (17) 



where p(R, Sl, Sr) is the probability density function, 
El(R, Sl, Sr) is the local energy, and Sl,r are the re- 
spective left (L), right (R) auxiliary, or shadow, variables. 
For the quartic wave function, the corresponding expres- 
sions for the probability density function and energy ex- 
pectation value are similar but with the addition of two 
more auxiliary shadow variables S' L R . 

We use the Metropolis Monte Carlo algorithm [21] 
to sample the probability density from a 9N and 15N- 
dimcnsional configuration space, corresponding to two 
and four sets of shadow coordinates in the case of lin- 
ear/quadratic and quartic wave functions respectively. 
In these computations, the Metropolis steps are subdi- 
vided in two parts. First, one attempts to move real 
particle coordinates at random inside cubical boxes of 
side length A. Second, analogous attempts are made 
to move shadow coordinates inside cubical boxes of side 
length A s h. For instance, in the case of the quadratic 
wave function, we first attempted to move all the R co- 
ordinates, then move the shadow coordinates {Sl}, then 
{Sr}. The parameters A and A s h were adjusted so that 
the acceptance ratio for both particle and shadow moves 
was nearly 50%. 

In addition to the ground state variational energy, 
we have also computed the radial distribution function 
g(r), and its Fourier transform, the structure factor S(k). 
These quantities are spherical averages and have been 
computed for both the real particles and the shadow co- 
ordinates. The radial distribution function is defined by 
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where the angular brackets denote an average with re- 
spect to \^f(R, t)| 2 and p is the particle density. The 
structure factor S(k) is obtained from the average 
jf(P-kPk), where p k = Ylf=i exp(-ik -rj), a procedure 
which is only possible on a discrete set of k values allowed 
by the periodic boundary conditions. 

All simulations presented in this work have been done 
with N — 108 atoms of 4 He in a cubic box with periodic 
boundary conditions. To enforce periodicity all correla- 
tions smoothly go to zero at a cutoff distance, r c = L/2, 



equal half the side of the simulation box according to the 
replacement 



/(r) - f(r) + f(2r c - r) - 2/(r c ) . 



(19) 
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FIG. 1. The Ground state energy per particle in Kelvin for 
4 He at the experimental equilibrium density (pa 3 = 0.365) 
using the Aziz HFDHE2 potential as a function of the pa- 
rameter t. Monte Carlo results from using various short-time 
evolved wave functions are as indicated. All simulations have 
been done for N = 108 particles. M indicates a McMillan 
wave function energy. M+MS, M+AS, M+T, OJ+AS refers 
to various variational Monte Carlo (VMC) results in the lit- 
erature, see text for details. 

In Fig. 1 we show the equilibrium energy per parti- 
cle for liquid 4 He for various short time evolved wave 
functions as function of the imaginary time parame- 
ter r. Other results from literature are also indicated 
for comparison: M + MS is the energy obtained by a 
shadow wave function having a pure repulsive McMil- 
lan (M) pseudopotential [22] of fifth power law form for 
both particles and shadows. [23] M+AS is the energy 
obtained by a shadow wave function with an attractive 
shadow-shadow pseudopotential of scaled Aziz HFDHE2 
potential (AS) form. [24] OJ+AS refers to a shadow 
wave function with an optimized Jastrow particle-particle 
pseudopotential (OJ) and scaled Aziz HFDHE2 shadow- 
shadow pseudopotential (AS). [24] GFMC is the Green's 
Function Monte Carlo calculations with Mcmillan form 
for importance and starting function. [4] The experimen- 
tal value is taken from Roach et al [25]. 

As expected, each of our factorized wave function 
reaches an energy minimum with increasing value of r. 
The flatness and depth of the energy minimun improve 
markedly with the order of the wave function. The lin- 
ear wave function has a shallow and narrow minimum 
at t = 0.002 and only improves upon McMillan's result 
(t = 0) by « 0.3 K. The minimum of the quadratic wave 
function is much better at r = 0.006 with a value of 
—6.393 K. The quartic wave function's energy minimum 



3 



extends further out to r = 0.015 attaining —6.809 K, 
which is lower than all existing variational Monte Carlo 
calculations that we are aware of. To demonstrate the ne- 
cessity of the gradient term, we have also plotted results 
obtained from (7) without the gradient term in the po- 
tential. In the present case, gradient term is responsible 
for w 50% of the improvement from that of the quadratic 
wave function. 

To give an quantatitive comparison, we summarize var- 
ious ground state equilibrium energies for 4 He in Table I. 

In Fig. 2 we show the equilibrium pair distribution 
function g(r) for 4 He as obtained from the quartic wave 
function. This g(r) is compared with the respective g(r) 
obtained from the M+AS shadow wave function and the 
experimental one of Svensson et al [26] obtained by neu- 
tron diffraction at saturated vapor pressure at T = 1.0 K. 
It is known [24] that the M+AS curve differs from the ex- 
perimental one because it predicts a diminished nearest- 
neighbor maximum and the entire curve is shifted by 
about 0.1 A to larger values of r compared to the exper- 
imental results. The pair distribution function that we 
obtained is in excellent agreement with the experimental 
one. 
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FIG. 2. The pair distribution function for liquid 4 He at 
the equilibrium density pa 3 = 0.365 after a VMC simulations 
with N=108 particles. The filled circles show the g(r) of this 
work that is compared with the respective g(r) obtained from 
the M+AS wave function (dotted line) and the experimental 
g(r) as reported by Svensson et al [26] (solid line) obtained 
at saturated vapor pressure at a temperature T — 1.0 K. 

In Fig. 3 we show S(k) at equilibrium density pa 3 — 
0.365 as obtained from the quartic wave function. The 
experimental S(k) shown in this figure is the result re- 
ported by Svensson et al [26]. The overall agreement 
between our short-time evolved structure factor with ex- 
periment is excellent except at small k. This is not un- 
expected because our imaginary time is still rather short 
for the wave function to develop the necessary long-range 



correlation to produce the linear behavior [27] of S(k) ob- 
served in bulk 4 He. 
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FIG. 3. Static structure factor S(k) of liquid 4 He at equilib- 
rium density pa s — 0.365. The filled circles show our results 
for S(k) obtained from the formula S(k) = ^-(p_kPk). The 
solid line denotes the experimental results reported by Svens- 
son and co-workers [26] obtained at saturated vapor pressure 
by means of neutron diffraction at temperature T = 1.0 K. 

In this work, based on recent findings on forward time 
steps decomposition schemes, we have implemented a 
fourth order Short-Time-Evolved wave function for de- 
scribing the ground state of strongly interacting quan- 
tum systems. Our approach is systematic, free of ar- 
bitrary parameters, and can be applied to any general 
quantum many-body problem. In the case of liquid 4 He, 
we have produced ground state energy and structure re- 
sults better than any existing VMC calculations, but 
without the use of complicated branching processes as 
in DMC or GFMC. Since the anti-symmetric require- 
ment on Fermion wave functions can be more easily im- 
plemented on the variational level, our quartic wave func- 
tions may be of great utility in studying Fermi systems. 

Our second order wave function is similar in structure 
to the class of shadow wave functions [28] , except that 
our wave function follows directly from the second order 
factorization scheme without any particular adjustment 
of pseudopotential or scale functions. Our use of a pos- 
itive factorization scheme to produce a much improved 
fourth order wave function demonstrates that there is 
a systematic way of improving this class of wave func- 
tions by introducting more shadow coordinates. Cur- 
rently, there is no known sixth order forward factorization 
schemes, and hence no sixth order Short-Time-Evolved 
wave function is possible. 
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TABLE I. Energies of liquid 4 He at the experimental equi- 
librium density (per 3 = 0.365 ; a = 2.556A) and at zero tem- 
perature. VMC indicates a variational Monte Carlo calcula- 
tion with the indicated wave function. All simulations use the 
Aziz HFDHE2 potential and have been performed for systems 
of N = 108 particles. The M+MS results are taken from Vi- 
tiello et al [23]. The M+AS and OJ-AS results are taken from 
MacFarland et al [24], The GFMC results are taken from Ka- 
los et al [4]. The experimental data are taken from Roach et 
al [25]. The energies are given in Kelvin per particle. 



Method 


Trial function 


Energy (K) 


VMC 


M+MS 


-6.061 ± 0.025 


VMC 


M+AS 


-6.599 ± 0.034 


VMC 


OJ+AS 


-6.789 ± 0.023 


VMC 


Linear 


-6.144 ± 0.092 


VMC 


Quadratic 


-6.393 ± 0.021 


VMC 


No Grad 


-6.644 ± 0.026 


VMC 


Quartic 


-6.809 ± 0.017 


GFMC 




-7.120 ± 0.024 


Experiment 




-7.140 
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